Steady-state mode I cracks in a viscoelastic triangular lattice 
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We continue our study of the exact solutions for steady-state cracks in ideally brittle viscoelastic 
lattice models by focusing on mode I in a triangular system. The issues we address include the 
' crack velocity versus driving curve as well as the onset of additional bond-breaking, signaling 

the emergence of complex spatio-temporal behavior. Somewhat surprisingly, the critical velocity 
for this transition becomes a decreasing function of the dissipation for sufficiently large values 
thereof. At the end, we briefly discuss the possible relevance of our findings for experiments on 
mode I crack instabilities. 
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^ I. INTRODUCTION 

I n the previous paper (Pechenik, ct al., 2000), we showed, how one could extend the methodology originally devised 
by Slepyan (1981 ) and obtain a closed form solution for mode III cracks propagating steadily in a square lattice in 
i_T the presence of dissipation in the form of a Kelvin viscosity. This approach allowed us to determine the effect of 
dissipation on the crack velocity as a function of driving. Similarly, we found the critical velocity (and the associated 
critical time) at which other bonds in the lattice have displacements larger than the assumed breaking criterion, 
signaling an instability in the steadily-propagating crack and the onset of more complex spatio-temporal behavior. 

In this second paper, we extend this analysis to the more interesting (and more experimentally relevant) case of 
mode I cracks. We now work on a triangular lattice and assume that there are central force ideally brittle vi scoelastic 



CNJ . springs connecting particles at the lattice sites. Again, we use ideas borrowed from Kulamekhtova (1984) to solve 
this problem exactly on an infinite lattice. We also solve the same system numerically on a lattice of finite transverse 
CD , extent so as to provide an independent check on some of the results. 

Perhaps the most interesting finding reported here concerns the aforementioned onset of additional bond-breaking 
at large enough velocity. At small dissipation, diagonal bonds that are offset from the assumed crack line in the vertical 
direction are the most "dangerous" . The critical velocity at which one of these bonds exceeds critical displacement is 
more or less independent of the dissipation at low values thereof and eventually rises sharply as the system becomes 
increasingly damped. This is similar to the behavior obtained for mode III. At moderate to large damping, however, 
horizontal bonds become more relevant and eventually give rise to a critical velocity which decreases with increasing 
damping. At the end, we will comment on the possible relevance of our results for trying to make sense of recent 



O ■ experiments ( Fineberg, ct al., 1991 , 1992 ; Bharon, ct al., 1995 , 199E) and simulations on instabilities during dynamic 
, fracture. 



II. MODE I ON A TRIANGULAR LATTICE 

Our model consists of central force springs between mass points located on a triangular lattice. We introduce unit 
vectors in the lattice directions as follows: d\ points in the direction of the x-axis and each of the subsequent unit 
vectors di with 2 < i < 6 rotated by an additional tt/3 in the counter-clockwise direction. This means that d i+3 = —di. 
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We will always assume that the steady-state crack breaks bonds between the rows with y = and y = v3/2, as it 
moves in the direction of the x-axis. 

As already discussed in our previous paper on mode III cracks, each bond is taken to be linear until a critical 
displacement at which point it completely breaks. Also, dissipation is incorporated as a Kelvin viscosity term. Given 
these assumptions, the equation describing the motion of the masses in the lattice is 

i— 1 

Here, x = (x,y) is a lattice vector in the plane and Qi(t,x) is the elongation of bond i emanating from the specific 
lattice site, in the approximation that this elongation is much smaller than the unstretched bond length, 

Qi(t, x) = (u(t, x + di) — u(t, x)j ■ di . (2) 

Finally, the driving term a(t, x) consists of forces which compensate for the forces from broken bonds (which are 
included on the left-hand side) as well as any external forces acting on the crack boundaries. 

Next we suppose that the crack is moving with a constant speed v. This allows us to change variables to r = x — vt, 
z = (r, y) , which gives 

^-Wg^^Qi^di-v 2 -^ = ~a(z) ■ (3) 
i=l 

Because of the form of Qi(z), we have the obvious symmetries 

Q 2 (i) = Q 5 {z + di), Q 3 {z) = Q e {z + d 3 ) . (4) 

We will impose one additional symmetry condition on the Qi. As the crack propagates, it alternately breaks bonds 
in the di and d 3 directions. Let us define our coordinates such that at r = and t = 0, (so that x = 0), the crack 
breaks a bond in the direction of d%. This means that for r > this bond (Qiizo) with zq = (r, 0)) is unbroken and 
for t < it is always broken. If we shift one lattice spacing to the left, the d^ bond at x = — 1 will break at t = — 1/v 
(here also t = in Q2{zq)). It is clear that the bond in the direction at x — breaks at a time-point between 
these two events; we will assume that this time is exactly at the midpoint, t = —l/(2v) or equivalently r = 1/2 in 
Qs(zq). This means that for r > 1/2 this bond is unbroken and for r < 1/2 it is always broken. We will assume an 
even more strict symmetry condition for these bonds, namely 

Q 2 (^o) = Qs(x - v(t - -L),y = 0) = Q 3 (z + £di) (5) 

which will need to be checked from our solution later. 

Let us denote Q%{zq) = Q(t). Then, we can write down an explicit expression for the forces B(z) on the right hand 
side of Eq. [| For the row with y — 

/ d 
a(z ) = 0(-t) I a e (z ) - (1 - Tjv—)Q 2 (z )d 2 

+0(-t + i) (a e '{3a) - (1 - V v^-)Q 3 (z )d^j (6) 

= *(-t) (-Po(r) - (1 - nv-jL)Q(r)) k 

+0(-r + \) (-P q (t - i) - (1 - riv-^Mr - ^ d 3 . (7) 

Here a e and a e ' are the assumed external forces. In the second form, these have been taken to match the vectorial 
nature and also the symmetry of the bond-breaking terms. The (scalar) Pq will be specified later. Note that the 
second theta-function in (||) reflects the already mentioned fact that the bond from point zq in the d 3 direction breaks 
earlier (by time interval l/2w) than the bond from the same point in the d 2 direction. 
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Let us write down the corresponding force for the row y = y/3/2. For the point x = di, the bond in the d§ direction 
is physically the same as the d 2 bond at x = (0,0) and therefore breaks at r = in Q§{zq + d 2 ) = Qziza)- Now, the 
de bond at x = d2 breaks later, at r = —1/2 in Qe(zQ + d®). Thus we have 

a(z + d 2 ) = 6(-t) (^ e (z + d 2 ) - (1 - Tqv^-)Q^(z + d 2 )d 5 ^j 

+0(-t - i) (a> e (z + d 2 ) - (1 - V v^-)Q 6 (z + d 2 )d^ (8) 

= 0(-t) (-P (t) - (1 - ^)Q(r)) d 5 + 

+0(~r \) (-Po(r + \) - (1 - Vv^)Q(r + i)) d 6 . (9) 

where we used the symmetry in Eq. || to write Qe(zo +d 2 ) = Qz(zq +d 2 — d%) — Q 2 (zq + d 2 — d% — d\/2) = Q(t + 1/2). 
Finally, we can rewrite all these forces in a more compact form, 

a(z ) = N(r)d 5 + N(t - i)<fe , (10) 

v{z a + d 2 ) = N{r)d 2 + N(t + ~)d 3 , (11) 

where 

N(r) = 9(-r) (p q (t) + (1 - r,v-^)Q(rfj . (12) 

Next, we proceed by doing a Fourier transformation of Eq. (^), over the continuous variable t and the discrete 
variable y — n^/Z/2, with integer n. In detail, 

u F (r,s) = g X,^nW S ", (13a) 



2 

71 — — CO \ 



u 



(T>V) = T- u F (T,s)e- ls yds; (13b) 



r+oo 

u FF (q,s) = I u F {T,s)e lqT dT , (14a) 



— oo 



We thereby obtain 



with 



1 f + °° 

u F (t,s) = — / u FF {q,s)e- lqT dq . (14b) 
2tt ./__, 



-v\^u FF - (1 + i V vq) J2 QFdi = ^ FF , (15) 



q ff = ( e -i?. di _ 1)u ff . £ (lg) 



and r = (q, s). 

To find a FF we first transform Eqs. ( |l0| , pi] ) over r 

^(g, ? i = 0) = iV(g)d 5 +e 4 37V(q)d 6 , (17a) 

3 F {q,n=l)e- l i = N(q)d 2 + N(q)d 3 , (17b) 

with 

N(q) = Pq + (1 + ir)vq)Q- - ifvQ(0) . (18) 
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The superscript "— " refers to the part of the Fourier transform which is anal ytic in the lower half plane of the variable 
q; our notation follows that in our previous mode III paper. Finally, using (13a) we get 



N{q)d 5 + e^N{q)d 6 + (e i iN(q)d 2 + N{q)d 3 



(19) 



This then completes the derivation of the basic equation that needs to be solved to determine the steady-state crack 
field. 



III. WIENER-HOPF SOLUTION 

We wish to use the Wiener-Hopf technique to solve the equation derived above for the elongation field Q. To start, 
we note that Eq. (15) can be written in matrix form 



M/ F = a FF 



(20) 



with the 2x2 matrix 



M 



A 1 A 3 
A 3 A 2 



(21) 



where 



Ai = -v 2 q 2 + (1 + irjvq) ^4 sin 2 

r ■ do 



A 2 = -v z q A + 3(1 + i-qvq) sin 



2 ' ' "2 . .2 

1 sin 



• 2 d 2 . 2 r ■ d 3 

sm h sm 

2 2 

r ■ d 3 \ 



A3 = V3(l + irjvq) sin 2 — sin 



2 r ■ d 2 . 2 f • d 3 



(22a) 
(22b) 
(22c) 



and 



N(q) ( (l-e 1 i)(l + e i ^ 5 ) 



y3(l + e^)(l-e l ^ s ) 



(23) 



The solution of (fed) is therefore 



^ FF _ 1 ( detMi 



detM V detM 2 



(24) 



where 



Mi 



a FF a FF 

w x y 

A 3 A 2 



Mo 



A x A 3 

FF FF 

U V y 



(25) 



Now, we want to extract from the solution (S3) an equation for the variable we are interested in, namely Q F (q). 
From the definition of Q we have 



Q F (q) = Q F (q,n = 0) 



^ ^(e^-l)^.4 ds 



4tt 



(26) 



From the basic solution above, we can evaluate u FF ■ d 2 



i FF -d 2 



d 2x A 3 + d 2y A x -d 2x A 2 + d 2y A 3 



.FF 



.FF 



IMI 



(27) 
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Before proceeding to evaluate the integral, it is worthwhile to pause and check the basic symmetry conditions. We 
find that 



»' ^ (e- ii! i *-l)u FF -d 3 ds, 



Q* (q,n = 0) = ^- 



(28) 



\/3 



where the new dot product is given by 



d 2x A 3 + d 2y A x d 2x A 2 + d 2y A 3 



T FF 



T FF 



IMI 



(29) 



Now, we can change variables of integration in (28 ) from s to — s. Under this transformation A\ 



A U A 2 



_ —iV3s/2 FF F F 



A 2 , A 3 \ s ^s = -A 3 , a FF \ s ^s 
and u FF ■ d 3 \ s ^- a = - e - l ^ s / 2 u FF . J 2 . Therefore ( ^8|) changes into 

V3 



°x F , Oy F \s^-s = -e- iVSs / 2 a FF , which gives that detM| s ^_ s = detM 



Q3 " = 0) 



4tt 



^qd^+isd 2y _ l) e ~id2yS U FF d2ds _ e HQ F ( q , n = 0) 



(30) 



This is exactly the symmetry condition (|5|); hence, our solution is consistent with the assumed symmetry. We also 
note here that a similar derivation gives 



1 1 v^3 

i(r,y) = Qz{t + -=,-y) = Qe(r + -,-y + — ) • 



(31) 



More generally, we can see what these symmetries mean for the displacement field u(z). Using detMij 



iV3*/2detMi, detM : 



8 s- 



= -e" lv/5s / 2 detM 2 , we obtain u FF \ 



= e-^ s l 2 u FF and u FF \ 



„ e -tV3s/2 u FF_ Thege giye u ^ y j = u ^ _ y + yg/ 2 ) y j = _ Uy (j^ _ y + ^3/2). 

We now proceed to calculate the integral in (|2^). The explicit form of the integrand in (|26| ) is 



N(q) 
IMI 



3 2f -(1 + irjvq){2 sin 2 § + sin 2 rfL ) „£JL + 3(1 + ir)vq ) sin 2 ^ 



(cosf • d% — 1) — (cos &s — cos |) (cosr • d 2 — 1) + ( cos -^r s — cos |) 



(32) 



As |M| is an even function of s, only the even part of the other determinant contributes to the integral. Then, by 
the change of variable tow = e -iV3s/2 we can rewr j^ e (26) as an integral over the unit circle in the complex w plane. 
After some tedious algebra, we find that we can write 



det M = 3£ + a£, + (3 



where 



£ = (l+irjvq)^ + w 

a = -2cos|(-2uV+3(l + ^wg)(l + 2sin 2 |) 
(3 = 3(1 + ir j vq f{l + 3 sin 2 |) - 4(1 + ir)vq)v 2 q 2 (l + sin 2 |) + v 4 q 4 

Similarly, the even part of the determinant in the numerator becomes 

F(0 



with 



no 



1 + irjvq ' 



(£ — (1 + iqrjv) cos ^) 2 + 2(1 + iqr]v) sin 2 -(1 + cos -)(1 + iqrjv — £) 



(33) 

(34a) 

(34b) 
(34c) 

(35) 



2 2 
-v q 



(1 + iqrjv — £)(1 + cos — ) — £ cos — + 1 + irjvq 



(36) 
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We write down here also the odd part of this determinant for future reference, 



G(w) = i( w) sin — (v 2 q 2 — 3(1 + iqrjv) sin 2 — ) 

w 2 2 



Inserting (B3J) and 



into ( |32j ) we can see that ( |26D becomes 
N(q) 



Q (q) 



no 



dw 



2wi(l + irjvq) J 3£ 2 + a£ + (3 w 



(37) 



(38) 



with integration around the unit circle in the counter-clockwise direction. The integrand in Eq. (|38|) has three poles 
inside the unit circle; i.e., w = and one w for each of the two roots of the expression (|33|), 



with 



a, 2 = c± Vc 2 - d 



C = cos^ I --v 2 q 2 + (1 + ir]vq)(l + 2sin 2 |) 



(1 + i?7ug) 2 (l + 3 sin 2 |) - ^(1 + irjvq)v 2 q 2 (l + sin 2 |) + i« 4 g 4 



Then calculating the residues at these poles we get 



Q^) = T ^L(l-5) 



vqvq 



where 



S(q) 



n&wtz - (! + ^ u «) 2 - nv*) vcf - (i + i^«) 2 



3(6 - 6)VS - (i + ^mPVli - (i + *w) 2 

Note that we are supposed to take those branches of the square roots in (^2|) which ensure that 

|wi,2| < 1 

with wi t 2 defined by 



Wl,2 



(1,2 - C 1 + *^9) 2 



irjvq 



(39) 

(40a) 
(40b) 

(41) 

(42) 

(43) 

(44) 



If we now substitute N(q) in Eq. ( 41 ) into Eq. (|l8|) , we get exactly the same equation as was obtained for the mode 
III case in the first paper (Pechenik, et al., 2000) 



Q+ + SQ- 



VV (1-S)Q(0)- {1 - S)P ° 



1 + ir/vq 



1 + irjvq 



(45) 



with Q having the same meaning as V in (Pechenik, ct al., 2000), namely the elongation of bonds between rows where 
the crack propagates. 

Therefore, we do not need to repeat any of the discussion from our previous paper regarding the choice of Pq" and 
the factorization into pieces analytic in the upper and lower half planes. Instead, we can just write down the final 
answer. For the relationship between the velocity and the driving displacement, we have 



where for this mode I case 



A _ / 1 + cjnyv 1 
A G ~ V A(f> 6XP 4iri 



A 



+ OC 




X(l + irrvx) 



(-l + v 2 f 



dx 



VsJl-v 2 



(46) 



(47) 
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FIG. 1 Dependence of A/A G on speed v for r] = 0.1, 0.7, 1.3, 1.9 
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FIG. 2 Change of behavior for small r\. Note the rapid oscillations at smaller speeds. 



and 



S 2 {q){q 2 + 4> 2 ) 



This expression for A contains the transverse speed Ct 



A 2 q 2 (f) 2 

= -\/3/8 and longitudinal speed ci 



KK1>= - "',"■>> ■ ( 48 ) 

/9/8. The zero of A 

gives the limiting speed of crack propagation c/j = 2 = 0-563, which is as expected equal to the Raleigh wave 
speed. 

In our previous paper, we discussed the details of how to evaluate this equation so as to obtain numerical results for 
the crack velocity response curve. These methods can be applied to this case as well and hence we do not repeat any 
of the details here. The results of these calculations are presented in Figs. [I] and ||. In fact, these findings are quite 
similar to the analogous mode III results. Later, we will show that the upper limiting value of A/ A<3 for an arrested 
(i.e. non-moving) crack is around 1.94. As was the case for mode III, the A — v curves in Figures [i] || approach this 
value as v goes to at any value of the damping r\. Slepyan's original calculations, done for the dissipationlcss limit 
r\ = 0, show the same asymptotic value. 



IV. CONSISTENCY OF THE STEADY-STATE SOLUTION 



In this section we calculate the bond displacements in the vicinity of the crack. This is crucial, as we need to check 
whether all bonds assumed to be unbroken in the derivation in fact have elongations less than unity, the value at 
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FIG. 3 Elongation of bonds on the crack surface, Q(r) for 77 = 0.5. 



which bonds break. 

Let us start with bonds along the crack path, where the elongation is just determined by transforming Q F (q) back 
to physical space. This part follows exactly the analogous mode III calculation, and we easily obtain for Q 



Q + = 
Q- = 



Q(o) 



s+ 



wQ(Q) 

(— iq + 0)(1 +iqr)v) 5+1 j_ 1 + iqrjv 

Q(Q) 1 | yvQ(0) 

S + \ „ = j_ S~ (iq + 0)(1 + iqrjv) 1 + iqrjv 



and for Q(t) 



Q(r) 
Q(r) 



Q(P) 



0(0) 



(K+)i dq 



— , for r > 



+00 



y/{iq - 0)(iq - 0)(1 + iqrfv) S+ \ q =^ ^ 

iq + <j)\^ (K^y^ dq 
2tt 



(49a) 
(49b) 

(50a) 



, y/A<j)(iq + Q)(1 + iqr)v) \iq + Q 
Q(0)e^, forr<0. 



'<}= — 



(50b) 



The numerical evaluation of these expressions follows the same methodology as described for mode III. Typical 
results are shown in Fig. [| We note that the elongation along the crack line is rather similar to the same object in 
the mode III case. In that case, however, this function also determined the horizontal bond elongations by simple 
subtraction. This is no longer true for mode I, since the vectorial nature of the problem requires that we take a 
different component of the displacement (different than the ones which goes into Q) to evaluate this elongation. In 
detail, we now have for 



Q h (q) = Q((q,n = 0) 



^1 (e- ir ^-l)u FF d 1 ds = 

47T J_ 2tl 

■S3 

V3 aetMi dg 
4tt J_2 JL detM 

V3 



(51) 
(52) 



where these matrices were defined in the last section. 

We now proc eed as before to change variables to w and re- write this expression in terms of the auxiliary variable £ 
defined in Eq. fl34a|) . 



detMi 



1 + i-qvq 



term odd in s 



(53) 
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FIG. 4 Elongation of the horizontal bonds, Qh(r), for r\ — 0.1 



where 



H{0 = f(l -e**)(l + ir)vq + i){\ + irjvq - ^ cos ^ - 
_| sin | ( l + e «f)(£2_ (1+^)2) 



(54) 



As before, the integrand in (52) has three poles. Dropping the irrelevant odd term and performing the integral leads 
to the expression 



Mq) = 



N(q) 
1 + ir/vq 



(1 



with 



i>h(q) = 



1 - e'i H{^)^e 2 - (1 + irrvW-H{h)yM - (1 + tyvq) 2 

2 3(6 - & - (i + ^) vel - (1 + «^) 2 



(55) 



(56) 



Branches of the square root satisfy the same conditions (^3|), as earlier. The function iV(q) can be found from 

Eq. @ 



(l + ^g)Q F (g) 
1-5 



Our final answer is 



Q(0) (1 



-iq 



S+\ 0= ^S~ (iq + 0)(l + iqr)v) 



(57) 



(58) 



For q close to 0, the function (1 — e~ lq )tjjh(q) behaves as q, giving a divergence 1/q 1 / 2 for Qh{q)- This is similar 
to the divergence in Q + (q). Thus, the numerical calculation of Qu(t) is similar to the calculation of Q(t) for r > 0. 
Figure ^ displays Qh(t) for several values of i] and v. Generally speaking, the function Qh(j) has maxima in two 
different places, one somewhere in the vicinity of —1 and a second for r > 0. 

We also need to find the bond elongation between the layers with n = 1 and n = 2 and between the layers with 
n = —1 and n = 0. Due to the symmetry (^lj) it is sufficient to consider only Q^iq^n = ±1), which we will denote 
as Q±(q). We can derive for it an expression similar to that in Eq. ( |38| ) with the major difference being that in this 
case the odd parts of the explicit determinant in Eq. (|32|) do not cancel out because of the additional factor w ±l . 



:(<?) 



±1 l+iyvq + G ( W ) dw 



N(q) _ 

2m f W 3f 2 + c< + /3 w 



N(q) I lS±GH 



2m J 3£ 2 + a£ + P 

\w\ = l 



dw 



(59) 
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FIG. 5 Elongation of bonds on the layer next to the crack layer, Q±(r) for rj = 0.1 



To derive the second equality, we used a change of variables w — ► 1/w, along with the symmetries of F and G. Again 
the integration in (^9|) is in the counter-clockwise direction. In this case, the integrand in the second form of the 
integral has just two poles inside the unit sphere, at W\ t 2. After the by now familiar tedious calculations, we find 



1 + irjvq 



(60) 



where 



with 



g(q)(w! - w 2 ) 

3(6 - 6) 



g(q) = -2i sin-(«V - 3(1 + iqrjv) sin 2 -) 



and N(q) is again given by (pffl. Finally we obtain 



Q±(q) 



Q(o) 



tp±(q) 



S+\ a= ^S- (iq + 0)(l+iqr]v) 



(61) 



(62) 



(63) 



Just as was the case for Qh(q) and Q + (q), Q±(q) behaves as 1/q 1 / 2 near 0. Thus numerical calculations can be 
performed just as in the previous cases. Figure |^ shows several curves Q±(t) for differing parameters. 

Given the elongations of these "vulnerable" bonds, we can investigate the critical speed at which one of these bonds 
should be broken. Figure || shows the results of our calculations of this critical speed. In Figure [7] we plotted r cr at 
which these bonds break; this will allow us to identify spatial and time coordinates of braking and make contact with 
our finite lattice calculations later. We found that the maximum value of Q-(t) always reaches 1 before Q+{t) and 
that this is the most dangerous bond for small dissipation. This curve turns around at an 77 of around 1.2, so that 
Q-(t) is always subcritical for larger rj. This is a result of the fact that the maximum value of Q_(r) is, surprisingly 
enough, not monotonic with velocity, but reaches a maximum and then decreases with increasing velocity. The 
maximum extension of this type of bond occurs, for large velocity, at some distance from the crack surface. For small 
rj, where the maximum Q_(t) does exceed critical extension, the decrease of the maximum Q_ with v restabilizes 
the bond beyond some velocity. The dominant threshold for rj above about 1.1, then, comes from the horizontal bond 
breaking. Note that as 77 increases, the critical velocity decreases. Also, for small i~j, a crossover occurs as the relative 
importance of the two different maxima in the horizontal bond elongation reverse; this region is anyway irrelevant as 
the next-layer vertical bond breaks at a lower velocity. Whenever the horizontal bond dominates, the maximum at 
t ~ — 1 is the governing one. As we see from Figure ^ even for large 77 this critical r stays near — 1, what means that 
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FIG. 7 Time of bond breaking at critical speed 



horizontal bond always breaks near the tip of the crack, contrary to the Mode III case, where the point of breaking 
drifts backward with increasing rj. 

Many of these features are special to the Mode I problem, and have no analog in the Mode III calculation. For Mode 
III, the horizontal bond breaking always dominates. Furthermore, the maximum off-crack surface bond extension is 
a strictly increasing function of the velocity. This points to the possibility that the dynamics of Mode I may be much 
richer than the Mode III dynamics. 



V. FINITE LATTICE MODEL. 



For our mode III calculations, it proved interesting to compare the exact results derived for an infinite lattice 



with the numerical deter mination of crack propagation properties in lattices of small transverse size (Kessler, 2000 
Kessler and Levine, 1998). We now discuss similar calculations for our mode I model. In addition to providing details 
regarding the size needed to attai n answers relevant to th e macroscopic limit (a question of direct relevance for direct 



molecular dynamics simu l ations (Abraham, et al., 1994 ; |Gumbsch, et al., 1997 ; Holland and Marder, 1997 , 1998 
Omeltchenko, et al., 1997; Zhou, et al, 1997), for example), finite lattice results can be used as a rough check on some 
of the predictions obtained above. Given the complexity of the analysis, having such a check is quite useful. 
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A. Arrested crack 

Let us start with trying to find arrested cracks, expected to exist in some range of the driving A around the 
Griffith's displacement. This involves looking for a solution of Eq. (|j) with u(x,y) independent of time and with 
the forces on the right hand side arising solely from the broken bonds. If we have a system with a finite number of 
rows, the natural boundary condition requires that the top and bottom rows have fixed vertical displacements of A 
and —A respectively. Since the entire system is linear, we can choose to use A = 1 and rescale the breaking criterion 
accordingly. 

Since we are doing a numerical calculation, we must introduce artificial boundaries in the direction along the crack, 
x. What we do is cut off the range of points whose coordinates are variables and outside of this range impose fixed 
asymptotic displacements. On the cracked side, these are just u = ± Ay for positive and negative y. For the uncracked 
side, the asymptotic displacement involves constant strain. At each of the sites that has a variable displacement, we 
impose the two components of the (vectorial) equation of motion. We also impose the equations of motion at the 
boundaries of the system. This gives us more equations than we have variables. The coupling of sites which contain 
variable displacements with the ones with fixed displacements gives rise to inhomogeneous terms. Combining these 
observations, we can write our system in the schematic form Mu — 6 = 0, with a non-square matrix M. The field 
u is then determined by the requirement that the error be minimal. In this manner, the small errors introduced at 
the boundary by having to have a fixed box size are prevented from causing any large errors (by modes which grow 
exponentially away from the edges) in the bulk of the lattice. 

The solution of this linear system determines the elongation of all bonds. In Figure ^ we show a typical lattice given 
by this solution. In fact we need to know elongation of just two bonds right at the tip of the crack; the one which 
in the case of the moving crack would break next (Q n ) and the other which would have been the last one broken 
(Qi). Now, recall that we have scaled our displacement to equal unity and the only remaining displacement scale is 
the elongation at which the springs break, which we can call e. The solution we have found is consistent as long as 
e is between the lower limit set by the Q n and the upper limit set by Qi. Since Aq = t^J (2N + l)/3, where 2N is 
the number of rows in the y-direction (excluding the boundary rows whose displacement is constrained), we directly 
obtain the upper and lower limits of the arrested crack band 

, 1 < A/A G < 1 (64) 

v /i(27V+l)Q ; v /i(27V+l)Q„ 

The results for these thresholds as a function of lattice size are shown in Figs. ^, [l0[ In the limit of infinite lattice, 
the lower limiting value A_/Aq approaches 0.515, while the upper limiting value A^/Aq approaches 1.94. 

B. Stable moving cracks 

We now turn to the moving crack problem. Again, we need to solve Eq. ([!]) with the forces in the right part due to 
the broken bonds. Now, the displacement field u(x, y, t) is of course time dependent. Therefore, we need to introduce 
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FIG. 9 Mode I on triangular lattice. Upper limiting A for arrested crack. Calculations were done for lattices with sizes 
2m x 2m. 
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FIG. 10 Mode I on triangular lattice. Lower limiting A for arrested crack. Calculations were done for lattices with sizes 
2m x 2m. 

a time-step At so as to define the time-points at which we will obtain a numerical solution. Given some specific speed 
v of the crack, we choose to divide half of the time interval that it takes for the crack to propagate one lattice spacing, 
l/(2v), into n equal time intervals; thus At — l/(2vn). Now we can discretize our system, obtaining equations at times 
t = 0, At, (n — \)At. We use a symmetric discretization of first and second derivatives (u(t + At) — u(t — At))/ (2 At) 
and (u(t + At) + u(t — At) — 2u(t))/ At 2 correspondingly. These equations depend on displacements u outside of this 
time interval, because of the temporal derivatives. These can be found via use of the assumed symmetries of the 
moving crack; in fact, it is easy to see that we thereby trade in displacements outside of the modeled interval with 
displacements inside this interval albeit at a different spatial location. The boundaries are treated the same way as 
described for the arrested crack; the displacements outside some range are replaced by asymptotic values for all of 
the time-points in our interval. Including the equations at these boundary points again gives us a system in which 
the number of equations exceeds the number of variables and again a least squared error algorithm is used to find the 
required solution. Fig. [ll] shows a snapshot of the lattice near the tip of the crack at one particular time for one set 
of parameters. 

We used these finite lattice calculations to provide an independent check on our analytic infinite lattice results. We 
were specifically interested in checking the critical speed estimates. The difficulty is that convergence in the parameter 
1/N is fairly slow; and, using large lattices is rather time-consuming and memory-demanding. In practice, we limited 
our calculations to TV's ranging from 6 to 24. It turns out that for small dissipation where the next-row vertical 
bond vulnerability is the thing which determines the critical velocity, we can do a credible job of verifying that the 
infinite lattice results are consistent with the finite lattice ones. For example, in Fig. we show the critical speed 
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FIG. 11 Snapshot of the region around tip of the crack for finite size model 2N x 2N for moving crack with v=0.48, r\ = 1, 
n = 5, N = 24. 
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as a function of lattice size for several different values of r). These numbers, if extrapolated to infinite N are clearly 
consistent with the results given earlier in Fig. ^ especially considering that having a finite time-step does introduce 
a small numerical error on its own. We can also check the critical time r at which the bond breaking occurs. From 
Fig. (0), we see that in the infinite lattice the break occurs for t ~ —0.5. For the finite lattice with n = 5, this mean 
that the bond a highlighted in Fig. [ll] should go above breaking threshold at the first time-point of the modeled 
interval. This is exactly what occurs. 

A different strategy is necessary for studying the onset of horizontal bond breaking. This is due to the fact that 
quite large A's are required in order to for the results to quantitatively approach the infinite N limit. This is another 
striking difference between the Mode I results and those of Mode III. For Mode III, the qualitative nature of the 
extraneous bond breaking was always similar to that of the infinite width limit, and the results were quantitatively 
accurate even for fairly small iV's. Here, however, the picture of which bonds are critical changes dramatically between 
finite N and the infinite N limit. In Fig. ([l3|), we present the "phase diagram" for critical bonds for N — 38, plotted 
together with the infinite N result, and in Fig. (|l4|), we compare the N = 38 results to those for N = 26. Since 
the steady-state code is difficult to run at these large values of N, the data in these two plots were obtained from 
time-dependent simulations, with only the central vertical bonds allowed to break. After running to times of t = 100 
to eliminate transients, the extensions of the dangerous bonds were checked to see if they exceeded criticality. Note 
that, in contradistinction to the "phase diagram" plotted in Fig. (||), we now plot A/Aq on the vertical axis, since 



15 




o.o 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 — 1 

0.0 0.5 1.0 1.5 2.0 2.5 3.0 



11 

FIG. 13 Comparison of the driving thresholds for critical bond extension for N — oo and N = 38 
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FIG. 14 Comparison of the driving thresholds for critical bond extension for N — 38 and iV = 26. 



this is the input control parameter for the simulations. We see that qualitative behavior of the bonds is similar 

to the infinite-TV result, with these bonds being below critical extension for both large rj and large A/Aq. The 
horizontal bonds also behave qualitatively like their infinite- N counterparts, but the quantitative agreement, is as we 
noted above, significantly worse. In fact, for large r/, the most dangerous bond is in fact of U Q+" type. However, 
the threshold at which the type bond becomes dominant is pushed to larger r\ as the system size is increased, 

presumably going to infinity with N. Also, there is a small region of 77 around 1.1 where the inconsistency is re-entrant, 
so that for intermediate values of A, no bonds of the crack surface are critical. Thus, to see dynamics qualitatively 
representative of the macroscopic system requires very large system sizes at moderate to large r\. Presumably, this 
is connected to the process zone increasing in size with r\. Again, it is worth emphasizing that this is a feature not 
present in Mode III, and leads to the conclusion that in molecular dynamics simulations, the width of the material 
should be taken very large to accurately study the micro-branching instability. 



VI. DISCUSSION 

In this work, we have discussed steady-state mode I crack propagation in a viscoelastic lattice model. Our primary 
method of analysis utilizes the Wiener-Hopf technique to write down closed form expressions for both the v — A 
curves (for various values of the dissipation 77) and for the bond elongation field. The latter enables us to define 
a limit of consistency for the solution past which some bonds not along the crack path have elongations greater 
than the assumed breaking criterion. These results are new, and correspond to non-trivial extensions of the work of 
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Kulamekhtova (1984) on the dissipationless limit and the work of Marder and Gross (1995) on small lattices. 

The most interesting results, in our opinion, concern the dependence of the critical velocity (where the aforemen- 
tioned inconsistency sets in) on the amount of dissipation. For small dissipation, this threshold is relatively insensitive 
to ?y, as already suggested in direct num erical simul ations. This threshold, occurring at roughly .73 of the Rayleigh 
speed is in some ways reminiscent of the Yoffe (1951) branching criterion which suggests that straight crack propaga- 
tion will become unstable once the largest stress direction shifts away from being straight ahead. It is unfortunately 
hard to be more precise regarding this correspondence since the crack dynamics in our model is fundamentally tied 
to the lattice scale, not the macroscopic scale - in fact, the latter is completely invisible in the Slepyan approach aside 
from providing a driving force in the form of a stress intensity factor. 

At larger 77, the instability picture changes. Now, it is a horizontal bond breaking which signals the onset of more 
complex crack dynamics. Also, the threshold goes down with increasing dissipation. This was not the case for our 
mode III calculations. This instability has nothing to do with Yoffe, as it is strongly dissipation dependent and in 
any case is not associated with crack branching. 

Much of the recent theoretical work on mode I cracks has been motivated by experiments which show clearly that 
instabilities limit the range of steady-state crack propagation. These instabilities introduce more complex spatio- 
temporal dynamics to the fracture process, causing additional dissipation and leaving behind a roughened crack 
surface. It has been tempting to associate these results with the onset of inconsistencies in lattice models, although 
the details of this correspondence remains uncertain. First, most of the experimental work has been carried out in 
amorphous materials, making the idea of an ordered lattice model somewhat suspect. Next, the instabilities seen 
experimentally typically occur at smaller speeds than the ones seen in lattice systems with small dissipation. Finally, 
the experiments seem to show a typical frequency for micro-branching which is not connected in any obvious way to 
dynamics at a lattice scale. 

Notwithstanding all these issues, we remain optimistic that the study of this class of models will lead to insight 
into dynamic fracture. There are several intriguing possibilities that need to be investigated in the future. First, we 
have shown that for mode I cracks, increasing the dissipation (in the form of a Kelvin viscosity) eventually results 
in a decrease of the instability threshold. If our model is really applied at the atomic scale, it is hard to see why 
there should be a large linear dissipation; on the other hand, it is well-known that lattice models miss an essential 
non-linear dissipation mechanism, namely the creation of dislocations which remain pinned to the crack. Would 
inclusion of these effects also lower the threshold? On the other hand, applying the model on a large scale (possibly 
for a disordered system) would naturally require a large dissipation and recent numerical simulations indicate that 



pro per inclusion of the thermal fluctua tions might also push the model into better agreement with experiment (Pla 
et a ., 199S ; Sander and Ghasias, 1999). Finally, the exact nature of the state which occurs past the instability onset 



has not been addressed in our work to date, and in fact cannot be addressed by the elegant but ultimately limiting 
analytic methods utilized for the steady-state problem. Instead, we plan to study a generalized force law in which the 
sharp breaking criterion is replaced by an analytic nonlinear force versus displacement ( Kessler and Levine, 1999 ). 
In t his formu lation, the inconsistency found here becomes a linear instability of the steady-state crack (Kessler and 
Lev ne, 2000 ) and one can use some of the methods developed in the field of nonequilibrium pattern formation to 
unravel the dynamics past onset. These studies, together with additional experimental data regarding the differences 
between brittle fracture in crystalline versus non-crystalline materials, will hopefully lead to a better understanding 
of dynamic fracture. 
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